
TString Y_Str[3];
int N_calltime = 0;
int NN_calltime = 0;

double m1 = 1115.68;    //MeV
// double m1 = 938.27;    //MeV
double m2 = 938.27;     //MeV
double m = m1*m2/(m1+m2);

double hbar = 197.3269631; // MeV*fm/c
double hbar_GeV = 0.1973269631; // GeV*fm/c

TFile* file_LL_momentum[3]; // Iy=0,1,2
TChain* chain_LL_momentum[3]; // Iy=0,1,2
TProfile* CF_LL[3]; // Iy=0,1,2
TVector3 *k[3]; // Iy=0,1,2
TVector3 *beta[3]; // Iy=0,1,2

double r_G[3][3];
double f_0,d_0;
TH1F* h_CF_map[3][3];
set<int> I_y_set;

double norm2(complex<double> f){
    double n = f.real()*f.real() + f.imag()*f.imag();
    // return sqrt(n);
    return n;
}

void GetCF_LL(TProfile* CF_LL_Fit,double rG,double f0,double d0,int I_y){ // pLambda
    CF_LL_Fit->Reset();
    double R_rel = rG*sqrt(2);
    Int_t Nevent = chain_LL_momentum[I_y]->GetEntries();
    // Int_t nevent = 250000;
    // cout<<"All Pair : "<<chain_LL_momentum[I_y]->GetEntries()<<endl;
    Int_t nevent = 125000;
    // Int_t nevent = 5000;
    // for(int n = 0;n<nevent;n++){
    for(int n = Nevent-1;n>=Nevent-nevent;n--){
        chain_LL_momentum[I_y]->GetEvent(n);
        TLorentzVector r_rel;
        r_rel.SetXYZT(gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel), 0);
        r_rel.Boost(*beta[I_y]);
        double r_star = r_rel.Vect().Mag(); // fm
        double k_star = k[I_y]->Mag()/hbar_GeV; // 1/fm
        double kr = k_star*r_star;
        double k_r = r_rel.Vect()*(*k[I_y])/hbar_GeV;
        double wei = 0;
        complex<double> C = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k_star*k_star,-k_star)); // fm
        wei = norm2( complex<double>(cos(k_r),-sin(k_r))+C*complex<double>(cos(kr)/r_star,sin(kr)/r_star) );
        CF_LL_Fit->Fill(k[I_y]->Mag()*1000, wei);
        // cout<<Form("k=%f,wei=%f\n",k[I_y]->Mag()*1000,wei);
    }
    return;
}


TH1F* h_CF_forFit;
void f0d0rG_fit_fcn(int& npar,double* deriv,double& f,double par[],int flag){
    NN_calltime++;
    // cout<<"\nNN_calltime="<<NN_calltime<<"\n"<<endl;

    double f0 = par[0];
    // double one_over_f0 = par[0];
    double d0 = par[1];
    double rG = par[2];
    double I_y_double = par[3];
    int I_y = int(round(I_y_double));

    CF_LL[I_y]->Reset();

    gRandom->SetSeed(4869); // 很重要！

    double R_rel = rG*sqrt(2);
    Int_t Nevent = chain_LL_momentum[I_y]->GetEntries();
    // Int_t nevent = 1000000 ;
    // Int_t nevent = 100000;
    // Int_t nevent = 250000;
    Int_t nevent = 125000;
    // Int_t nevent = 25000;
    // Int_t nevent = 20000;
    // Int_t nevent = 10000;
    // Int_t nevent = 5000;
    // Int_t nevent = 1000;
    // for(int n = 0;n<nevent;n++){
    for(int n = Nevent-1;n>=Nevent-nevent;n--){
        chain_LL_momentum[I_y]->GetEvent(n);
        TLorentzVector r_rel;
        r_rel.SetXYZT(gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel),gRandom->Gaus(0,R_rel), 0);
        r_rel.Boost(*beta[I_y]);
        double r_star = r_rel.Vect().Mag(); // fm
        double k_star = k[I_y]->Mag()/hbar_GeV; // 1/fm
        double kr = k_star*r_star;
        double k_r = r_rel.Vect()*(*k[I_y])/hbar_GeV;
        double wei = 0;
        complex<double> C = complex<double>(1,0)/(complex<double>(1/f0+0.5*d0*k_star*k_star,-k_star)); // fm
        wei = norm2( complex<double>(cos(k_r),-sin(k_r))+C*complex<double>(cos(kr)/r_star,sin(kr)/r_star) );
        CF_LL[I_y]->Fill(k[I_y]->Mag()*1000, wei);
    }
    int N_bin = h_CF_forFit->GetNbinsX();
    f = 0;
    // for(int n_bin = 1;n_bin<=N_bin;n_bin++){
    for(int n_bin = 1;n_bin<=20;n_bin++){
        double y_data = h_CF_forFit->GetBinContent(n_bin);
        double y_data_err = h_CF_forFit->GetBinError(n_bin);
        double y = CF_LL[I_y]->GetBinContent(n_bin);
        double y_err = CF_LL[I_y]->GetBinError(n_bin);
        // f = f + (y-y_data)*(y-y_data)/(y_data_err*y_data_err + y_err*y_err);
        f = f + (y-y_data)*(y-y_data)/(y_data_err*y_data_err);
        // cout<<Form("n_bin=%d,f=%f,y-y_data=%f,y_data_err=%f\n",n_bin,f,y-y_data,y_data_err);
    }
    cout<<Form("NN_calltime=%d,f=%f\t",NN_calltime,f);
    cout<<Form("\t\t f0=%.6f,d0=%.6f,rG=%.6f\n",f0,d0,rG);


    return;
}


double f_min,fedm,errdef;
int npari,nparx,istat;
double all_chi2_LL_TM(double f_0,double d_0){
    // NN_calltime++;
    N_calltime++;
    cout<<"\nNN_calltime="<<NN_calltime<<"\n"<<endl;
    cout<<"\nN_calltime="<<N_calltime<<"\n"<<endl;

    double chi2 = 0;
    
    for(int cent_i = 1; cent_i<4;cent_i++){  //// loop cut_ii <<<
        for(int I_y=0;I_y<3;I_y++){
            if(I_y_set.find(I_y) == I_y_set.end()){continue;}
            
            h_CF_forFit = h_CF_map[cent_i-1][I_y];

            const int npar = 4;
            TMinuit minuit(npar);
            minuit.SetPrintLevel(-1);//-1,0,1
            minuit.SetFCN(f0d0rG_fit_fcn);
            // minuit.DefineParameter(0,"f0",f_0,0.1,0,3.5);
            minuit.DefineParameter(0,"f0",f_0,0.1,0,4);
            // minuit.DefineParameter(1,"d0",d_0,0.001,-9.5,9.5);
            // minuit.DefineParameter(1,"d0",d_0,0.1,-9.5,9.5);
            minuit.DefineParameter(1,"d0",d_0,0.1,-6,6);
            minuit.DefineParameter(2,"r_G",3,0.1,0.01,5);
            minuit.DefineParameter(3,"I_y",I_y,0.5,0.01,5);
            minuit.FixParameter(0);
            minuit.FixParameter(1);
            minuit.FixParameter(3);

            minuit.SetErrorDef(1);
            // minuit.SetErrorDef(9);
            minuit.SetMaxIterations(500);
            minuit.Migrad();
            double rG,rG_err;
            minuit.GetParameter(2,rG,rG_err);
            
            minuit.mnstat(f_min,fedm,errdef,npari,nparx,istat);

            r_G[cent_i-1][I_y] = rG;
            
            chi2 = chi2 + f_min;

            cout<<Form("r_G[%d][%d] = %f , f_min=%f , chi2_now=%f\n",cent_i-1,I_y,r_G[cent_i-1][I_y],f_min,chi2);
        }
    }
    cout<<Form("f0=%f,d0=%f,all_chi2=%f\n,",f_0,d_0,chi2);
    return chi2;
}

void f0d0_fit_fcn(int& npar,double* deriv,double& f,double par[],int flag){
    double f0 = par[0];
    // double one_over_f0 = par[0];
    double d0 = par[1];
    // f = all_chi2(f0,d0);
    // f = all_chi2_LLNew(f0,d0);
    f = all_chi2_LL_TM(f0,d0);
    // f = all_chi2(1/one_over_f0,d0);
    return;
}

void Get_f0d0_TMinuit(double &x1,double &y1,double &x1_err,double &y1_err){
    // const int npar = 3;
    const int npar = 2;
    TMinuit minuit(npar);
    minuit.SetPrintLevel(-1);//-1,0,1
    // minuit.SetPrintLevel(0);//-1,0,1
    minuit.SetFCN(f0d0_fit_fcn);

    minuit.DefineParameter(0,"f0",2,0.1,0,4.5);
    minuit.DefineParameter(1,"d0",3.11026e+00,0.1,-6,6);

    minuit.SetErrorDef(1);
    // minuit.SetErrorDef(4);
    // minuit.SetErrorDef(9);
    minuit.SetMaxIterations(500);
    // minuit.SetPrintLevel(-1);//-1,0,1
    minuit.Migrad();
    minuit.GetParameter(0,x1,x1_err);
    minuit.GetParameter(1,y1,y1_err);

    cout<<Form("f0_min=%f+-%f,d0_min=%f+-%f,chi2_min=%f\n",x1,x1_err,y1,y1_err,all_chi2_LL_TM(x1,y1));


    return;

    auto file = new TFile("IyAll_f0_d0_20230220.phase3.5k.2.root","recreate"); 
    // auto file = new TFile("IyAll_f0_d0_20230215.1w.root","recreate"); 
    TCanvas* c2 = new TCanvas("Canvas_IyAll_f0_d0","Canvas_IyAll_f0_d0",1000,1000);

    gPad->SetMargin(0.12,0.12,0.12,0.12);
    gPad->SetFrameLineWidth(3);
    gStyle->SetLineWidth(3);

    double Error2 = 1;
    // auto contour_name=Form("IyAll_1overf0_d0_Contour%.2f",Error2);
    auto contour_name=Form("IyAll_f0_d0_Contour%.0f",Error2);
    minuit.SetErrorDef(Error2);
    gStyle->SetTitleFontSize(0.09);
    TGraph *gr = (TGraph*)minuit.Contour(100,0,1);
    // TGraph *gr = (TGraph*)minuit.Contour(10,0,1);
    // gr->SetTitle("Chi2 Contour of d_{0} and f_{0}  (-0.5<y<0);f_{0}/fm;d_{0}/fm");
    gr->SetTitle("Chi2 Contour of d_{0} and f_{0}  (-1<y<0);f_{0}/fm;d_{0}/fm");
    gr->SetName(contour_name);
    gr->SetFillColor(42);
    // gr->Draw("alf");
    gr->SetMinimum(-8);
    gr->SetMaximum(10);
    // gr->SetAxisRange(1.5,3);

    gr->GetXaxis()->SetLimits(0.5,4.5);

    gr->GetXaxis()->SetTitleSize(0.05);
    gr->GetXaxis()->SetTitleOffset(0.9);
    gr->GetXaxis()->SetLabelSize(0.035);
    gr->GetYaxis()->SetTitleSize(0.05);
    gr->GetYaxis()->SetTitleOffset(0.7);
    gr->GetYaxis()->SetLabelSize(0.035);

    gr->SetLineWidth(3);
    gr->Draw("al");
    gr->Write();

    gr->Print();
    c2->Print("test20230220_IyAll.1.phase3.5k.2.png");
    c2->Print("test20230220_IyAll.1.phase3.5k.2.root");


    Error2 = 4;
    // contour_name=Form("IyAll_1overf0_d0_Contour%.2f",Error2);
    contour_name=Form("IyAll_f0_d0_Contour%.0f",Error2);
    minuit.SetErrorDef(Error2);
    gr = (TGraph*)minuit.Contour(100,0,1);
    // TGraph *gr = (TGraph*)gMinuit->Contour(100,0,1);
    gr->SetName(contour_name);
    gr->SetFillColor(42);
    // gr->Draw("alf");
    gr->SetLineWidth(3);
    gr->Draw("same");
    gr->Write();

    gr->Print();
    c2->Print("test20230220_IyAll.2.phase3.5k.2.png");
    c2->Print("test20230220_IyAll.2.phase3.5k.2.root");

    Error2 = 9;
    // contour_name=Form("IyAll_1overf0_d0_Contour%.2f",Error2);
    contour_name=Form("IyAll_f0_d0_Contour%.0f",Error2);
    minuit.SetErrorDef(Error2);
    gr = (TGraph*)minuit.Contour(100,0,1);
    // TGraph *gr = (TGraph*)gMinuit->Contour(100,0,1);
    gr->SetName(contour_name);
    gr->SetFillColor(42);
    // gr->Draw("alf");
    gr->SetLineWidth(3);
    gr->Draw("same");
    gr->Write();

    gr->Print();
    c2->Print("test20230220_IyAll.3.phase3.5k.2.png");
    c2->Print("test20230220_IyAll.3.phase3.5k.2.root");

    TGraphErrors *gre_marker = new TGraphErrors(1);
    gre_marker->SetName("IyAll_f0_d0");
    gre_marker->SetMarkerColor(kRed);
    gre_marker->SetLineColor(kRed);
    gre_marker->SetPoint(0,x1,y1); gre_marker->SetPointError(0,x1_err,y1_err);
    gre_marker->SetLineWidth(3);
    gre_marker->Draw("pe");
    gre_marker->Write();

    c2->Print(TString("IyAll_f0_d0_20230220.phase3.5k.2")+TString(".png"));
    c2->Write();

    

    file->Close();

    return;
}



void LL_Fit_20230215(){

    Y_Str[0] = TString("(-0.5<y<0)");
    Y_Str[1] = TString("(-0.75<y<-0.5)");
    Y_Str[2] = TString("(-1<y<-0.75)");


    TString Cent_String[4];

    Cent_String[0] = TString("(Cent 0-60%)");
    Cent_String[1] = TString("(Cent 0-10%)");
    Cent_String[2] = TString("(Cent 10-20%)");
    Cent_String[3] = TString("(Cent 20-60%)");

    for(int I_y =0;I_y<3;I_y++){
        chain_LL_momentum[I_y] = new TChain(Form("pLambda_Momentum_Iy%d",I_y));
        chain_LL_momentum[I_y]->Add(Form("pLambda_Momentum_Iy%d.5kPair1Bin.3.root",I_y)); // 6000w
    
        k[I_y] = new TVector3(1,2,3);
        beta[I_y] = new TVector3(1,2,3);
        chain_LL_momentum[I_y]->SetBranchAddress("k", &(k[I_y]) );
        chain_LL_momentum[I_y]->SetBranchAddress("beta", &(beta[I_y]) );//beta
        Int_t nevent = chain_LL_momentum[I_y]->GetEntries();
        cout<<Form("nevent for Iy%d=%d\n",I_y,nevent);
        
        CF_LL[I_y] = new TProfile(Form("CF_pLambda_Iy%d",I_y), Form("CF_pLambda_Iy%d",I_y), 25, 0, 250);
    }

    // auto file_data = new TFile("/data4/qinz18/STAR/KF_LambdaProton/Result/0503/LambdaProton_AllCF_20220323_0503_Final_CF_SmearCorrection_HistBook.root","READ");
    auto file_data = new TFile("LambdaProton_AllCF_20220323_0503_Final_CF_SmearCorrection_HistBook.root","READ");


    for(int cent_i = 1; cent_i<4;cent_i++){  //// loop cut_ii <<<
        for(int I_y=0;I_y<3;I_y++){
            TH1F* hLambdaProton_CF = h_CF_map[cent_i-1][I_y];
            h_CF_map[cent_i-1][I_y] = (TH1F*)((TH1F*)file_data->GetObjectChecked(Form("h_Lambda_CF_Cent%d_CUT%d",cent_i,I_y), "TH1F"))->Clone();
            hLambdaProton_CF = h_CF_map[cent_i-1][I_y];
            hLambdaProton_CF->GetXaxis()->SetTitle("k/MeV");hLambdaProton_CF->GetYaxis()->SetTitle("C(k)");
            hLambdaProton_CF->SetMinimum(0.4);
            hLambdaProton_CF->SetMaximum(3.6);
            hLambdaProton_CF->SetAxisRange(0,200.5);
        }
    }

    double f0_TM;
    double d0_TM;
    double f0_TM_err;
    double d0_TM_err;

    gStyle->SetLineWidth(3);
    gStyle->SetTitleFontSize(0.09);

    I_y_set = {0,1,2};
    // I_y_set = {0};

    if(0){// Get TH2D for chi2 <<<
        int nbin = 20;
        double f0_min = 0;
        double f0_max = 4.5;
        double d0_min = -5;
        double d0_max = 5;

        TH2D* h2_chi2 = new TH2D("h2_chi2","Chi2 for pLambda LL_Fit;f0/fm;d0/fm",nbin,f0_min,f0_max,nbin,d0_min,d0_max);
        for(int i1=0;i1<20;i1++){
            double f0_i = f0_min+i1*(f0_max-f0_min)/20;
            for(int i2=0;i2<20;i2++){
                double d0_i = d0_min+i2*(d0_max-d0_min)/20;
                double chi2 = 0;
                for(int I_y=0;I_y<3;I_y++){
                    for(int cent_i = 1; cent_i<4;cent_i++){
                        h_CF_forFit = h_CF_map[cent_i-1][I_y];

                        const int npar = 4;
                        TMinuit minuit(npar);
                        minuit.SetPrintLevel(-1);//-1,0,1
                        minuit.SetFCN(f0d0rG_fit_fcn);
                        minuit.DefineParameter(0,"f0",f0_i,0.1,0,4.5);
                        minuit.DefineParameter(1,"d0",d0_i,0.1,-6,6);
                        minuit.DefineParameter(2,"r_G",3,0.1,0.1,5);
                        minuit.DefineParameter(3,"I_y",I_y,0.5,0.1,5);
                        minuit.FixParameter(0);
                        minuit.FixParameter(1);
                        minuit.FixParameter(3);

                        minuit.SetErrorDef(1);
                        // minuit.SetErrorDef(9);
                        minuit.SetMaxIterations(500);
                        minuit.Migrad();
                        double rG,rG_err;
                        minuit.GetParameter(2,rG,rG_err);
                        
                        minuit.mnstat(f_min,fedm,errdef,npari,nparx,istat);
                        chi2 = chi2 + f_min;
                    }
                }
                int bin_i = h2_chi2->FindBin(f0_i,d0_i);
                h2_chi2->SetBinContent(bin_i,chi2);
            }
        }

        TCanvas* c_chi2 = new TCanvas("c_chi2","c_chi2",2200,1200);
        h2_chi2->Draw("colz");
        c_chi2->Print("pLambda_h2_chi2.png");
        c_chi2->Print("pLambda_h2_chi2.root");

    }// Get TH2D for chi2 >>>

    // return;


    Get_f0d0_TMinuit(f0_TM,d0_TM,f0_TM_err,d0_TM_err);
    cout<<Form("f0_TM=%.10f,d0_TM=%.10f,f0_TM_err=%.10f,d0_TM_err=%.10f\n",f0_TM,d0_TM,f0_TM_err,d0_TM_err);
    // return;

    TCanvas* c1 = new TCanvas("c1","c1",2200,1200);
    TPad* drawpad = new TPad("drawpad", "drawpad",0.01,0.01,0.99,0.95);

    gPad->SetMargin(0.12,0.12,0.12,0.12);
    // gPad->SetMargin(0.2,0.2,0.2,0.2);
    drawpad->SetMargin(0.2,0.2,0.2,0.2);
    // gPad->SetFrameLineWidth(3);
    gPad->SetFrameLineWidth(4);

    drawpad->SetGrid();
    // drawpad->Divide(4,3,0.001,0.001);
    // drawpad->Divide(3,3,0.001,0.001);
    drawpad->Divide(3,3,0,0);
    drawpad->Draw();
    // drawpad->SetFillStyle(1);
    drawpad->SetFillStyle(0);
    // drawpad->SetTitle(TString("RawCF"));
    // drawpad->SetTitle(TString("CF_pLambda_UrQMD_Crab_LLFit"));
    drawpad->SetTitle(TString("CF_pLambda_Data_LLFit"));

    for(int I_y=0;I_y<3;I_y++){
        double chi2 = 0;
        // I_y_set = {I_y};
        I_y_set = {0,1,2};
        // Get_f0d0_TMinuit(f0_TM,d0_TM,f0_TM_err,d0_TM_err);
        // cout<<Form("f0_TM=%f,d0_TM=%f,f0_TM_err=%f,d0_TM_err=%f\n",f0_TM,d0_TM,f0_TM_err,d0_TM_err);
        for(int cent_i = 1; cent_i<4;cent_i++){  //// loop cut_ii <<<
            drawpad->cd(3*I_y+cent_i);

            h_CF_map[cent_i-1][I_y]->Draw();
            // h_CF_map[cent_i-1][I_y]->SetTitle((Cent_String[cent_i]+Y_Str[I_y]).ReplaceAll(")("," , ").ReplaceAll("("," ").ReplaceAll(")"," "));
            h_CF_map[cent_i-1][I_y]->SetTitle("");


            h_CF_map[cent_i-1][I_y]->GetXaxis()->SetTitleSize(0.08);
            h_CF_map[cent_i-1][I_y]->GetXaxis()->SetTitleOffset(0.9);
            h_CF_map[cent_i-1][I_y]->GetXaxis()->SetLabelSize(0.05);
            h_CF_map[cent_i-1][I_y]->GetYaxis()->SetTitleSize(0.09);
            h_CF_map[cent_i-1][I_y]->GetYaxis()->SetTitleOffset(0.5);
            h_CF_map[cent_i-1][I_y]->GetYaxis()->SetLabelSize(0.07);

            if(I_y!=2){
                h_CF_map[cent_i-1][I_y]->GetXaxis()->SetTickLength(0);
                // h_CF_map[cent_i-1][I_y]->GetXaxis()->SetLabelOffset(999);
            }
            if(cent_i!=1){
                h_CF_map[cent_i-1][I_y]->GetYaxis()->SetTickLength(0);
                // h_CF_map[cent_i-1][I_y]->GetYaxis()->SetLabelOffset(999);
            }

            h_CF_forFit = h_CF_map[cent_i-1][I_y];

            const int npar = 4;
            TMinuit minuit(npar);
            minuit.SetPrintLevel(-1);//-1,0,1
            minuit.SetFCN(f0d0rG_fit_fcn);
            minuit.DefineParameter(0,"f0",f0_TM,0.1,0,4);
            minuit.DefineParameter(1,"d0",d0_TM,0.1,-6,6);
            minuit.DefineParameter(2,"r_G",3,0.1,0.01,5);
            minuit.DefineParameter(3,"I_y",I_y,0.5,0.01,5);
            minuit.FixParameter(0);
            minuit.FixParameter(1);
            minuit.FixParameter(3);

            minuit.SetErrorDef(1);
            // minuit.SetErrorDef(9);
            minuit.SetMaxIterations(500);
            minuit.Migrad();
            double rG,rG_err;
            minuit.GetParameter(2,rG,rG_err);
            
            minuit.mnstat(f_min,fedm,errdef,npari,nparx,istat);

            r_G[cent_i-1][I_y] = rG;
            
            chi2 = chi2 + f_min;

            cout<<Form("r_G=%f,f0=%f,d0=%f,chi2=%f,cent%d,y%d\n",rG,f0_TM,d0_TM,f_min,cent_i,I_y);

            TProfile* CF_LL_Fit = new TProfile(Form("CF_pLambda_Fit_Cent%d_Iy%d",cent_i,I_y), Form("CF_pLambda_Fit_Cent%d_Iy%d",cent_i,I_y), 25, 0, 250);
            GetCF_LL(CF_LL_Fit,rG,f0_TM,d0_TM,I_y);
            CF_LL_Fit->SetLineColor(kRed);
            CF_LL_Fit->Draw("same");


            TPaveText* pt_title = new TPaveText(0,3.25,200,3.6,"");//中上
            pt_title->AddText((Cent_String[cent_i]+Y_Str[I_y]).ReplaceAll(")("," , ").ReplaceAll("("," ").ReplaceAll(")"," "));
            pt_title->SetFillStyle(0);
            pt_title->SetBorderSize(0);
            pt_title->Draw("same");


            TPaveText* pt = new TPaveText(100,1.8,180,3,"");//右上
            pt->AddText(Form("r_{G} = %.3f #pm %.3f",rG,rG_err));
            pt->AddText(Form("f_{0} = %.3f #pm %.3f",f0_TM,f0_TM_err));
            pt->AddText(Form("d_{0} = %.3f #pm %.3f",d0_TM,d0_TM_err));
            pt->Draw("same");

        }
        cout<<Form("total chi2 = %f for %s\n",chi2,Y_Str[I_y].Data());
    }

    c1->cd();
    TPad *padtitle = new TPad("padtitle", "padtitle",0.01,0.95,0.99,1);
    padtitle->Draw();
    padtitle->cd();
    padtitle->SetFillStyle(0);

    TString title("p#Lambda Correlation Function (LL Fit)");

    TPaveText* pt = new TPaveText(0,0,1,1);
    pt->AddText(title);
    pt->SetFillStyle(0);
    pt->SetBorderSize(0);
    pt->Draw();

    c1->Print("pLambda_f0_d0_All_20230215.png");
    c1->Print("pLambda_f0_d0_All_20230215.root");
}